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Magnetic Order and Charge Disproportionation in a Spin-ice type Kondo Lattice Model: Large 

Scale Monte Carlo Study 
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Phase diagram of a spin-ice type Kondo lattice model, potentially relevant to metallic pyrochlore oxides, 
is obtained by the Monte Carlo simulation implementing the polynomial expansion technique up to the 
system size with 2048 sites. We identified a new 32-sublattice magnetic phase with concomitant charge 
disproportionation, along with other phases such as two-in two-out and all-in/all-out orders. The spin 
and charge pattern can be switched by external magnetic field to a different one accompanied by a half 
magnetization plateau. 
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Geometrical frustration offers a fertile ground for study- 
ing fascinating phenomena in strongly correlated systems. 1,2) 
Competing interactions under the frustration often lead to 
extensive number of energetically-degenerate states. Even a 
small perturbation to the degeneracy can result in dramatic 
effects, such as phase transitions and colossal responses to 
external fields. This sensitivity has stimulated intensive stud- 
ies of competing orders and fluctuations in the geometrically- 
frustrated systems. 

A representative example of frustrated systems is the 
spin ice. 3_5) In these systems, spins with strong Ising-type 
anisotropy along the sublattice-dependent local (111) direc- 
tion reside on the pyrochlore lattice, which consists of corner- 
sharing tetrahedra, as shown in Figs. l(c)-(f). The nearest- 
neighbor (NN) interaction between the Ising spins is domi- 
nantly ferromagnetic (FM), favoring the two-in two-out con- 
figuration in each tetrahedron. This two-in two-out constraint, 
called the ice rule, leads to macroscopic degeneracy and pre- 
vents the system from ordering. 6) The situation offers an inter- 
esting playground for controlling the peculiar magnetic states 
by external magnetic field. 5 ' 7) In addition, it was theoretically 
predicted that the magnetic states are drastically changed 
by the sign and strength of exchange interactions. 8) Namely, 
when the dominant interaction changes from FM to antiferro- 
magnetic (AFM), all-in or all-out spin configuration becomes 
favored in each tetrahedron, resulting in a long-range order of 
alternating all-in and all-out tetrahedra. Although such phase 
transition is intriguing, the sign of interaction is usually inher- 
ent to each material and not flexibly tunable. 

In the present study, we explore another route to tuning the 
ice-type frustrated systems. We here consider a coupling of 
the spin-ice type moments to itinerant electrons. The kinetic 
motion of electrons causes an effective interaction between 
the localized moments, called the Ruderman-Kittel-Kasuya- 
Yosida (RKKY) interaction. 9) The effective interaction is, in 
general, long-ranged and oscillating (can be both FM and 
AFM), sensitively depends on the electron density and band 
structure. Therefore, the system has great potential to tune the 
phase competition by controlling carrier doping as well as ex- 
ternal pressure which modifies the band structure. 



Potential realization of such situation is brought by metallic 
pyrochlore oxides, such as 7?2Mo207 and R2lv20i (R is rare- 
earth), in which Ising-like rare-earth / moments couple with 
itinerant d electrons. The interplay between spin and charge is 
believed to play an essential role in various peculiar features 
in these systems, such as the anomalous Hall effect 10 ' n) and 
resistivity minimum. 12) 

In this Letter, we elucidate thermodynamic properties of an 
Ising- spin Kondo lattice model on the pyrochlore lattice by 
employing an unbiased numerical simulation. A related study 
was performed for a spin-only model with the long-range 
RKKY interaction. 13) However, the RKKY form is based on 
a second-order perturbation in terms of the spin-charge cou- 
pling with assuming a simple isotropic Fermi surface, neglect- 
ing the actual band structure. In the present study, we solve the 
problem in the nonperturbative region by a large scale Monte 
Carlo (MC) simulation which includes the electronic degree 
of freedom explicitly, and reveal a new spin-charge coupled 
phase and its switching by external magnetic field. 

We consider a Kondo lattice model with Ising spins on a 
pyrochlore lattice, whose Hamiltonian is given by 



(i,j),o- 



(1) 



The first term is hopping of itinerant electrons, where (cj^) 
is the annihilation (creation) operator of an itinerant electron 
with spin cr =f, I at zth site, and t is the transfer integral. The 
sum (i, j) is taken over NN sites on the pyrochlore lattice. The 
second term is the onsite interaction between localized spins 
and itinerant electrons, where S; and &{ represent the local- 
ized Ising spin and itinerant electron spin at zth site, respec- 
tively (| Sf | = 1), and / is the coupling constant (the sign of 
/ does not matter in the present model as the localized spins 
are classical). The anisotropy axis of Ising spin is given along 
the local (111) direction at each site, i.e., along the line con- 
necting the centers of two tetrahedra which the spin belongs 
to. Hereafter, we take t - 1 as the unit of energy, the lattice 
constant of cubic unit cell a = 1 [see Fig. 1(f)], and the Boltz- 

mann constant = 1. 

The model in eq. (1) was studied by the authors by using a 
cluster dynamical mean-field theory with a focus on transport 
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properties. 14) In the previous study, they calculated the phase 
diagram with limiting the magnetic unit cell to a single tetra- 
hedron. It is highly desired to make a complete phase diagram 
by considering longer-range spatial correlations. 

Here we study thermodynamic properties of the model by 
an unbiased MC simulation. The partition function of the 



present model is given by Z = Tr^jTr 



,exp[-{tf({S;})- 



fiN c }/T] (N c is the fermion number operator and ji is the 
chemical potential). The former trace is taken over the Ising 
spin degree of freedom, which is handled by the MC sam- 
pling using the standard Metropolis method with single spin 
flip dynamics. To obtain the MC weight, the latter trace for 
itinerant electrons is calculated by the following two meth- 
ods. One is the exact diagonalization (ED), 15) which is appli- 
cable only to small system sizes because the computational 
amount increases rapidly as a function of the system size 
N as 0(N 3 ). The other method is the polynomial expansion 
method (PEM), 16) which reduces the computational amount 
to 0(N 2 logA/), 17) enabling us to access larger system sizes. 
Prior to calculations, we tested the convergence of PEM in 
terms of the order of Chebyshev polynomials m. 18) In rel- 
atively higher density region n e > 0.15, we concluded that 
m = 40 gives good convergence (n e = ^io-( c Jo- c io-)/N). On the 
other hand, for n e < 0. 15, the convergence appears to be much 
worse, and hence we employed ED. Most of the calculations 
were conducted up to N = 4 x 6 3 by PEM and iV = 4x4 3 
by ED for typically more than 3000 MC steps. In addition, 
calculations for N - 4 x 8 3 were also done at a filling in the 
32-sublattice ordered region. Note that this is a state-of-the- 
art calculation; one MC step for N = 4 x 8 3 takes about 50 
seconds by using 1024 CPU cores in the System B (SGI Altix 
ICE 8400EX) at ISSP supercomputer center. 

Figure 1(a) shows the phase diagram of the model (1) 
at / = 2 obtained by the MC simulation. There are four 
dominant phases at low T in the calculated density region 
< n e < 0.3. The open symbols show the critical tem- 
peratures T c for each phase determined from the inflection 
point or sharp jump of T dependence of the order parameter 
at each system size (see Fig. 2 below). The order parameter 
is defined by M q = [S^ max /N t ] l/2 , where S^ max is the maxi- 
mum component of sublattice-dependent spin structure factor, 



Jq* = ZijeaiSi • Sj) exp [iq • (r; - Tj)] /N t (N t = N/4 is the 
number of tetrahedra, a denotes four sublattices in a primi- 
tive cell, r ; is the position of ith site, and q is the wave num- 
ber). The inflection point gives a good estimate of T c , since, 
in general, the order parameter monotonically goes to zero (a 
nonzero value) as increasing N in the higher (lower) T region; 
the system size dependence of T c is reasonably small, while 
we do not have enough sets of data to perform the systematic 
finite-size scaling analysis. 19) The ordering pattern below T c 
is determined from the characteristic wave number q together 
with the local spin correlations within tetrahedra introduced 
below. 

In the low density region n e < 0.04, the system develops the 
q = two-in two-out order below T c < 0.025 [Fig. 1(c)]. We 
call this phase the ice-ferro phase hereafter. While increasing 
the electron density, a different ordering develops in the region 
of 0.08 <n e < 0.15. 20) The ordering structure is of layer type 
with q = (0, 0, 27r); every tetrahedron retains two-in two-out 
ice-rule configuration, while the net moments of tetrahedra 
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Fig. 1. (color online), (a) Phase diagram of the model (1) at / = 2. The 
open symbols show T c for the four phases, while the closed ones are for 
other complex orders. 20 ^ The lines are guides for the eyes. The bottom 
stripe shows the ground-state phase diagram obtained by comparing the 
energy of four phases. 2T> PS indicates phase separation, (b) RKKY inter- 
actions for the nearest-, second-, and third-neighbor spins calculated by us- 
ing the second-order perturbation in terms of / for the model (1). The latter 
two are multiplied by factor 3 for clarity, (c)-(f) show spin configurations 
for (c) ice-ferro, (d) ice-(0, 0, 2n), (e) 32-sublattice, and (f) all-in/all-out 
order, respectively. 



form a collinear layer-type ordering [Fig. 1(d)]. We call this 
the ice-(0, 0, 2n) phase. The ordering pattern is the same as 
that found in a spin ice model including the long-range dipole- 
dipole interaction. 8) On the other hand, in the high density 
region n e > 0.22, the system exhibits the q = all-in/all-out 
order [Fig. 1(f)]. 

In the intermediate region of 0.15 < n e < 0.22, another 
complex magnetic structure appears. The magnetic ordering 
has a 32-sublattice structure characterized by the wave vector 
q = (7r,7r,7r). 21) The complicated structure can be viewed as 
a collection of one-dimensional chains with "in-in-out-out" 
spin configuration [see Fig. 1(e)]; namely, all the second- 
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neighbor spins along the chains are AFM. Surprisingly, such 
chain AFM is free from frustration in the pyrochlore lat- 
tice, uniquely selecting the 32-sublattice ordering. In terms 
of tetrahedra, the 32-sublattice order consists of a periodic ar- 
rangement of all 2 4 = 16 possible tetrahedra, as shown in 
Fig. 1(e): six two-in two-out, four three-in one-out, four one- 
in three-out, one all-in, and one all-out type. In this compli- 
cated structure, the mean-field from the neighboring spins is 
spatially inhomogeneous; namely, the mean-field from NN 
spins cancels to zero only at the sites belonging to all-in 
and all-out tetrahedra, while that at other sites has the same 
nonzero magnitude. This point will be discussed later in re- 
lation to charge disproportionation and its switching by mag- 
netic field. 

We also calculated the phase diagram at T = by a varia- 
tional calculation with limiting the trial states to the four or- 
dered states appearing in the MC simulation. 22) The result is 
shown at the bottom axis of Fig. 1(a). All the four phases ap- 
pear in the corresponding density regions where the MC result 
shows their instabilities. The phase transitions between differ- 
ent magnetic phases are all first order at T = 0, and hence, 
they are accompanied by a jump of n e , that is an electronic 
phase separation (PS). 23) 

Let us discuss the phase diagram from the viewpoint of the 
RKKY interactions. Figure 1(b) shows the nearest-, second-, 
and third-neighbor components of the RKKY interaction cal- 
culated by using the eigenstates for the / = tight-binding 
model. In the lowest density region n e < 0.05, all three com- 
ponents are positive, namely FM, because of the small Fermi 
surface. This is consistent with the ice-ferro order [Fig. 1(c)]. 
By increasing the electron density above 0.05, the second- and 
third-neighbor RKKY interactions change their signs to be 
AFM, but the NN interaction remains dominantly FM. The 
ice-(0, 0, 2ti) order is stabilized as a compromise of these in- 
teractions [Fig. 1(d)]. With further increase of the electron 
density, the NN RKKY interaction also changes its sign at 
n e ~ 0.20. The 32-sublattice order emerges in this region 
where the NN interaction becomes irrelevant; indeed, it is 
characterized as AFM ordering of the second-neighbor spins 
along the chains, as mentioned above. In the higher density 
region n e > 0.22, the NN interaction becomes AFM, and the 
all-in/all-out order is stabilized. 

It should be stressed that, although the sequence of mag- 
netic orderings is qualitatively understood, it is still hard to 
tell the phase diagram only from the RKKY interactions. In 
particular, the complicated 32-sublattice order is difficult to be 
predicted by the simple RKKY analysis. 13) The critical tem- 
perature of each phase is also difficult to be predicted; for in- 
stance, although the magnitude of NN ferromagnetic RKKY 
interaction at n e ~ 0.1 is comparable to that of the antifer- 
romagnetic one at n e ~ 0.3, T c is much lower in the former 
region compared to the latter. 14) The unbiased MC simulation 
explicitly taking account of itinerant electron degree of free- 
dom is crucial to identify the phase diagram and T c . 

Next, let us look into T dependences of the order parame- 
ters and short-range spin correlations within tetrahedra. The 
local spin correlations are measured by calculating the ratio 
of all-in or all-out (P40), three-in one-out or one-in three-out 
(P31), and two-in two-out (P22) tetrahedra. Figure 2(a) shows 
the results in the ice-ferro region. The results show an en- 
hancement of the ice-rule local correlation P22 prior to phase 
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Fig. 2. (color online). T dependences of the order parameter M q and the 
ratios of tetrahedra with different spin configurations, P22, P31, an d P40, 
at J = 2. The data are calculated at a fixed chemical potential ju: (a) fi = 
-5.9 [corresponding to n e = 0.030(2)], (b) p = -4.8 [n e = 0.099(3)], (c) 
fi = -3.4 [n e = 0.195(1)], and (d) fi = -2.4 [n e = 0.264(3)]. The magnetic 
ordering wave vectors for each phase are (a),(d) q = 0, (b) (0, 0, 2n), and 
(c) (n,n,7t). In (c), the charge disproportionation n q with q = (0, 0, 2n) is 
also plotted. 



transition at T c ^ 0.023, while P40 and P31 are strongly sup- 
pressed. P22 becomes dominant also in the ice-(0, 0, 2n) phase 
region, as shown in Fig. 2(b). 24) On the other hand, the results 
for all-in/all-out region show contrastive behavior [Fig. 2(d)], 
in which the all-in/all-out correlation P 40 is enhanced prior to 
the phase transition at T c 0.051. In sharp contrast, for the 
32-sublattice order, the local correlation parameters show lit- 
tle T dependences, even in the critical region near T c ^ 0.047. 
This implies that the ordering is not driven by NN spin corre- 
lations, as discussed above. 

Furthermore, the 32-sublattice magnetic order exhibits a 
concomitant charge disproportionation. In Fig. 2(c), we plot 
the charge disproportionation defined by n q = [N^ m3LX /N t ] 1/2 . 

Here, N ( q a) = Z iJea <w,-n/> ex P [iq ' ( r * - / N t is the charge 
structure factor and rti - YjA^i^io) i s me l° ca l electron den- 
sity at each site. The result clearly indicates the emergence of 
concomitant charge disproportionation below T c (the temper- 
ature of the inflection point of n q agrees with that of M q ). The 
wave number is characterized by q = (0, 0, 2n); the local elec- 
tron density is higher at the sites belonging to the all-in/all-out 
tetrahedra compared to the other sites. The pattern appears to 
be related with the inhomogeneity of the mean- fields from 
NN spins mentioned above. We note that all the other phases 
in Fig. 1 are charge uniform, and the 32-sublattice order is 
the only phase showing a concomitant charge disproportiona- 
tion. 2 ^ 

The peculiar charge disproportionation can be switched by 
external magnetic field. Figure 3 shows the magnetization 
M = I Yji Sfl/iV and charge disproportionation under the mag- 
netic field h applied along the [111] direction at T = 0.025. 
Here, the charge disproportionation An is defined by the dif- 
ference of the electron densities between the kagome and tri- 
angular planes perpendicular to the [111] direction. For sim- 
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also electronic and transport properties induced by the spin- 
charge interplay in frustrated itinerant electron systems. 
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Fig. 3. (color online). Magnetization M and charge disproportionation An 
under the external magnetic field h along the [111] direction at \i = -3.4 
and T = 0.025. The vertical dashed line is an estimate of the critical field 
at T = 0, h c 0.0465, obtained by comparing the ground state energy 
in a sufficiently large system. The inset shows a schematic picture of the 
magnetic order in the plateau state at M 1/4. 

plicity, we apply the magnetic field only to the localized Ising 
spins. Under weak magnetic field up to h ~ 0.04, the 32- 
sublattice order at h = remains robust and M stays almost 
zero. For higher magnetic field, M increases abruptly, indicat- 
ing a first order transition to a different phase, which is char- 
acterized by the half magnetization plateau M - 1/4. (The 
full saturation is M = 1/2 in the present model.) This plateau 
state remains stable up to h ~ 0.1. The magnetic structure of 
the plateau phase is obtained by aligning all the spins on the 
triangular layers to the field direction (see the inset of Fig. 3). 

Interestingly, An also changes from almost zero to nonzero 
abruptly at the critical field, namely, the charge dispropor- 
tionation is simultaneously switched to the one along the 
[111] direction. 26) The switching may be understood as fol- 
lows. As mentioned above, the charge disproportionation with 
q = (0, 0, 2n) in the zero-field state appears to be dominantly 
driven by the inhomogeneity of the NN spin correlations. In 
contrast, in the plateau state in applied magnetic field, the NN 
mean-fields have uniform, nonzero magnitude at all the sites. 
Instead, the mean-fields from second-neighbor spins along the 
chains are different between the triangular and kagome sites, 
which may lead to the [1 1 1] charge disproportionation in the 
plateau state. 

To summarize, we have numerically investigated the Kondo 
lattice model with Ising spins with (111) anisotropy on a 
pyrochlore lattice, which potentially describes metallic py- 
rochlore oxides. By the state-of-the-art Monte Carlo simu- 
lation with taking account of the itinerant electron degree 
of freedom, we found the 32-sublattice ordered phase in the 
competing region between the two-in two-out and all-in/all- 
out phases. This phase exhibits a charge disproportionation 
concomitant with the magnetic order. The spin and charge pat- 
tern can be switched by external magnetic field to the differ- 
ent one which shows a half-magnetization plateau. Our result 
demonstrates that the spin-charge coupling on the frustrated 
lattice induces richer behaviors than in the localized spin sys- 
tems. It will stimulate further study of not only magnetic but 
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24) The sharp jump of M q in Fig. 2(b) is indicative of a first order transition. 
This is presumably due to the six-fold degeneracy of the ordered ground 
state. In this respect, the transition to the ice-ferro state might also be of 
first order, while it is not clear in Fig. 2(a) presumably because of the 
small system size. 

25) A standard second-order perturbation in terms of J predicts that the un- 
derlying magnetic ordering restricts a possible wave number for charge 
disproportionation to q = q ai + q ai + G, where q a are the magnetic 
wave numbers and G is a reciprocal lattice vector. For the ice-ferro, 



stripe, and all-in/all-out phases, no charge disproportionation is allowed 
because 2q a = G. For the 32-sublattice order, however, charge dispro- 
portionation with q + G is allowed as q a is dependent on a. 
26) Although there is a finite size effect, the data for N = 4 x 6 3 show 
the transition to the plateau state very close to h c estimated at T = 
in the thermodynamic limit. In the high field region above the plateau 
state, however, we could not obtain converged results due to the poor 
convergence of PEM. 



